Seasonal dietary changes influence gut microbiota composition depending on the host species but do not alter gut microbiota diversity in lizards
Mauricio Hernández
Laboratory of Biotic Interactions, Centro de Investigación en Ciencias Biológicas, Universidad Autónoma de Tlaxcala, Tlaxcala, Mexicodelmer.hernandez@unah.edu.hn
Sergio Ancona
Instituto de Ecología, Universidad Nacional Autónoma de México, Mexicosergio.ancona@iecologia.unam.mx
Stephanie Hereira-Pacheco
Laboratory of Biotic Interactions, Centro de Investigación en Ciencias Biológicas, Universidad Autónoma de Tlaxcala, Tlaxcala, MexicoAníbal H. Díaz de la Vega-Pérez
Universidad Autónoma de Tlaxcala, Tlaxcala, MexicoAntton Alberdi
Center for Evolutionary Hologenomics, The GLOBE Institute, University of Copenhagen, Copenhagen, DenmarkYendi E. Navarro-Noya
Laboratory of Biotic Interactions, Centro de Investigación en Ciencias Biológicas, Universidad Autónoma de Tlaxcala, Tlaxcala, Mexicoyendiebenezer.navarro.n@uatx.mx
29 - 03 - 2024
I. PLOTS OF TAXONOMIC COMPOSITION
Loading libraries
PLOT RELATIVE ABUNDANCE AT THE PHYLUM LEVEL
Load data
Create phyloseq objects
Data visualizations
## [1] "Sg.F.D.1" "Sg.F.D.2" "Sg.M.D.3" "Sg.M.D.4" "Sg.F.D.5" "Sg.F.D.6"
## [7] "Sg.F.D.7" "Sg.F.D.8" "Sg.F.D.9" "Sg.M.D.10" "Sg.F.D.11" "Sg.M.D.12"
## [13] "Sg.M.D.13" "Sg.M.D.14" "Sg.M.D.15" "Ss.M.D.16" "Ss.M.D.17" "Ss.M.D.18"
## [19] "Ss.F.D.19" "Ss.M.D.20" "Ss.F.D.21" "Ss.M.D.22" "Ss.M.D.23" "Ss.M.D.24"
## [25] "Ss.F.D.25" "Ss.M.D.26" "Sa.M.D.27" "Sa.F.D.28" "Sa.M.D.29" "Sa.F.D.30"
## [31] "Sa.M.D.31" "Sa.M.D.32" "Sa.F.D.33" "Sa.F.D.35" "Sa.M.D.36" "Sb.F.D.37"
## [37] "Sb.M.D.38" "Sb.F.D.39" "Sb.M.D.41" "Sb.M.D.42" "Sb.M.D.43" "Sb.F.D.44"
## [43] "Sb.F.D.45" "Sb.M.D.46" "Ss.M.R.52" "Ss.F.R.53" "Ss.M.R.54" "Ss.M.R.55"
## [49] "Ss.F.R.56" "Ss.F.R.57" "Ss.F.R.58" "Ss.M.R.59" "Ss.M.R.60" "Ss.M.R.61"
## [55] "Sg.M.R.62" "Sg.M.R.63" "Sg.M.R.64" "Sg.M.R.65" "Sg.M.R.66" "Sg.F.R.67"
## [61] "Sg.F.R.68" "Sg.F.R.69" "Sg.M.R.70" "Sg.M.R.71" "Sg.F.R.72" "Sg.M.R.74"
## [67] "Sg.M.R.75" "Sa.F.R.72" "Sa.F.R.73" "Sa.M.R.74" "Sa.M.R.75" "Sa.F.R.76"
## [73] "Sa.M.R.77" "Sa.M.R.78" "Sa.M.R.79" "Sa.F.R.80" "Sa.M.R.81" "Sa.F.R.82"
## [79] "Sa.F.R.83" "Sa.F.R.84" "Sb.M.R.85" "Sb.M.R.86" "Sb.M.R.87" "Sb.F.R.88"
## [85] "Sb.F.R.89" "Sb.F.R.90" "Sb.M.R.91" "Sb.F.R.92" "Sb.M.R.93" "Sb.M.R.94"
## [91] "Sb.M.R.95" "Sb.F.R.96" "Sb.F.R.97" "Sb.M.R.98" "Sb.F.R.99"## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"## [1] "Season" "Species_Scel" "Sex"
## [4] "SeqDepth" "Weight" "SVL"
## [7] "TL" "Tb" "Ta"
## [10] "Ts" "Elevation" "Location"
## [13] "Date" "BarcodeSequence" "LinkerPrimerSequence"RELATIVE ABUNDACE BY TAXONOMIC LEVEL (MEAN AND STANDARD DEVIATION)
PLOT RELATIVE ABUNDANCE AT THE PHYLUM LEVEL
Create color palette
## [1] "#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C" "#FB9A99" "#E31A1C" "#FDBF6F"
## [8] "#FF7F00" "#CAB2D6" "#6A3D9A" "#FFFF99" "#B15928" "gray" "orange"Species + Season
Rel_Abun_Phyl_Season <- plot_bar(physeq = relative, "Sample", fill = "Phylum") +
facet_grid(~Species_Scel + Season, space = "free", scales = "free") +
labs(y="Relative abundance") +
geom_bar(stat = "identity", pisition = "stack", res = 300) +
scale_fill_manual(values = paleta) +
theme(strip.text.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold")) +
theme(strip.text.x = element_text(face = "italic")) +
theme(legend.text = element_text(face = "italic")) +
theme(strip.background = element_rect(colour = "black", fill = "white")) +
theme(text = element_text(size = 10))## Warning in geom_bar(stat = "identity", pisition = "stack", res = 300): Ignoring
## unknown parameters: `pisition` and `res`
PLOT RELATIVE ABUNDANCE AT THE GENUS LEVEL
Load data
metadata <- read.csv(file = "../data/metadata.csv", check.names = F)
asv_table <- read.csv("../data/feature_table.csv", check.names = F)
taxonomy <- read.csv("../data/taxonomy_03.csv", check.names = F) %>% unite(
taxa, Kingdom:Species, remove = F, sep = ";")
asvtable_metadata <- asv_table %>%
inner_join(taxonomy)## Joining with `by = join_by(OTUID)`Relative abundance at phylum level (Sceloporus species)
Genus_01 <- asvtable_metadata %>%
group_by(Genus) %>% # But also for other taxonomic levels
summarise_if(is.numeric, sum)
Genus_01 <- Genus_01 %>%
column_to_rownames(var = "Genus")
Genus.ra <- t(t(Genus_01)/colSums(Genus_01)*100)
rowMeans(Genus.ra) %>%
as.data.frame() %>% arrange(desc(.))## Alistipes Bacteroides
## 2.748686 8.221501
## Clostridia_vadinBB60_group Hafnia-Obesumbacterium
## 2.075043 8.502137
## Hungatella Lachnoclostridium
## 3.843127 1.646522
## Odoribacter Others
## 3.661410 12.190726
## Parabacteroides Pseudomonas
## 3.747141 5.199502
## Roseburia [Eubacterium]_coprostanoligenes_group
## 3.872776 1.770494Genus <- Genus.ra %>% t() %>%
as.data.frame()%>%
rownames_to_column(var = "SampleID") %>%
inner_join(metadata)## Joining with `by = join_by(SampleID)`Mean_Genus <- Genus %>%
group_by(Species_Scel) %>%
summarise_if(is.numeric, mean)
SD_Genus <- Genus %>%
group_by(Species_Scel) %>%
summarise_if(is.numeric, sd)
aggregate(Genus[ ,2:13], list(Genus$Species_Scel), mean)Filter genus taxonomic level
metadata <- read.csv(file = "../data/metadata.csv", header = TRUE, row.names = 1)
asv_table <- read.csv(file = "../data/feature_table.csv", check.names = F)
taxonomy <- read.csv("../data/taxonomy_03.csv", check.names = F) %>% mutate_at(
c("Genus"), str_replace,"g__", "")lista <- rowMeans(Genus.ra) %>% as.data.frame() %>%
arrange(desc(.)) %>%
slice_head(n=14) %>%
rownames_to_column(var = "Genus") %>%
filter(!Genus =="g__") %>%
mutate_at(c("Genus"), str_replace,"g__", "")
list <- lista$Genus
lista01 <- read.csv(file = "../data/lista.csv", check.names = F)
list02 <- lista01$Genus
#write.table(lista, file="../data/lista.txt", sep = "\t")## Joining with `by = join_by(Genus)`taxonomy_1 <- taxonomy_filter %>%
inner_join(asv_table, by =c(
"OTUID"="OTUID")) %>% dplyr::select(1:8)Load data again
Create phyloseq objects
Convert to relative abundance
Color palette
## [1] "#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C" "#FB9A99" "#E31A1C" "#FDBF6F"
## [8] "#FF7F00" "#CAB2D6" "#6A3D9A" "#FFFF99" "#B15928" "#1B9E77" "#D95F02"
## [15] "#7570B3" "#E7298A" "#66A61E" "#E6AB02" "#A6761D" "#666666"Plot by Species and Season
Rel_Abun_Genus_Season <- plot_bar(physeq = relative, "Sample", fill = "Genus") +
facet_grid(~Species_Scel + Season, space = "free", scales = "free") +
labs(y="Relative abundance") +
geom_bar(stat = "identity", pisition = "stack", res = 300) +
scale_fill_manual(values = paleta) +
theme(strip.text.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold")) +
theme(strip.text.x = element_text(face = "italic")) +
theme(legend.text = element_text(face = "italic")) +
theme(strip.background = element_rect(colour = "black", fill = "white")) +
theme(text = element_text(size = 10))## Warning in geom_bar(stat = "identity", pisition = "stack", res = 300): Ignoring
## unknown parameters: `pisition` and `res`
Core microbiota (50%)
aeneus_50 <- read.delim("../data/core_50_aeneus.tsv", sep = "\t", check.names = F) %>%
rownames_to_column(var = "ids")
bicanthalis_50 <- read.delim("../data/core_50_bicanthalis.tsv", sep = "\t", check.names = F) %>%
rownames_to_column(var = "ids")
grammicus_50 <- read.delim("../data/core_50_grammicus.tsv", sep = "\t", check.names = F) %>%
rownames_to_column(var = "ids")
spinosus_50 <- read.delim("../data/core_50_spinosus.tsv", sep = "\t", check.names = F) %>%
rownames_to_column(var = "ids") Plotting the Venn Diagramm
venn.plot_50 <- venn.diagram(
x = list(S_grammicus = grammicus_50$OTUID,
S_spinosus = spinosus_50$OTUID,
S_aeneus = aeneus_50$OTUID,
S_bicanthalis = bicanthalis_50$OTUID),
category.names = c(
expression(bold("S. grammicus")),
expression(bold("S. spinosus")),
expression(bold("S. aeneus")),
expression(bold("S. bicanthalis"))),
filename = "../figures/viendo_50.jpg",
output = TRUE,
height = 3000,
width = 6000,
resolution = 300,
compression = "lzw",
units = "px",
lwd = 6,
lty = "blank",
fill = c("#56B4E9", "#009E73", "#999999", "#E69F00"),
cex = 1.5,
#fontface = "bold",
fontface = "italic",
fontfamily = "sans",
cat.cex = 2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.pos = c(-10, 25, -0.020, -0.050),
cat.dist = c(0.3, 0.35, 0.099, 0.099),
cat.fontfamily = "sans")II. ALPHA DIVERSITY ANALYSES
Loading libraries
Calculating Hill numbers with hilldiv R package
Preparing the data
Load files
Declare factors
Evaluate whether sex variable influence on gut microbiota diversity (q=1)
# Paired test (Wilcoxon test, q1)
# Sceloporus aeneus
aeneus_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" &
Sex == "Female", q1, drop = TRUE)
aeneus_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" &
Sex == "Male", q1, drop = TRUE)
aeneus_FemMale_GutMic <- wilcox.test(x= aeneus_fem_GutMic, y= aeneus_male_GutMic)
aeneus_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: aeneus_fem_GutMic and aeneus_male_GutMic
## W = 49, p-value = 0.4779
## alternative hypothesis: true location shift is not equal to 0# Sceloporus bicanthalis
bica_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" &
Sex == "Female", q1, drop = TRUE)
bica_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" &
Sex == "Male", q1, drop = TRUE)
bica_FemMale_GutMic <- wilcox.test(x= bica_fem_GutMic, y= bica_male_GutMic)
bica_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: bica_fem_GutMic and bica_male_GutMic
## W = 69, p-value = 0.9095
## alternative hypothesis: true location shift is not equal to 0# Sceloporus grammicus
gram_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" &
Sex == "Female", q1, drop = TRUE)
gram_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" &
Sex == "Male", q1, drop = TRUE)
gram_FemMale_GutMic <- wilcox.test(x= gram_fem_GutMic, y= gram_male_GutMic)
gram_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: gram_fem_GutMic and gram_male_GutMic
## W = 109, p-value = 0.5676
## alternative hypothesis: true location shift is not equal to 0# Sceloporus spinosus
spi_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" &
Sex == "Female", q1, drop = TRUE)
spi_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" &
Sex == "Male", q1, drop = TRUE)
spi_FemMale_GutMic <- wilcox.test(x= spi_fem_GutMic, y= spi_male_GutMic)
spi_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: spi_fem_GutMic and spi_male_GutMic
## W = 55, p-value = 0.8596
## alternative hypothesis: true location shift is not equal to 0Differences in snout vent length (SVL) among lizard species
set.seed(1234)
Data_SVL <- subset(Micro_div, select = c(Species_Scel, SVL))
kruskal.t <- Data_SVL %>% kruskal_test(SVL ~ Species_Scel)
kruskal.tVisualization: box plots with p-values
Post_hoc <- Post_hoc %>%
add_xy_position(x = "Species_Scel")
ggboxplot(Data_SVL, x = "Species_Scel", y = "SVL", fill = "Species_Scel") +
xlab(element_blank())+
scale_fill_manual(values = c("#56B4E9","#009E73", "#999999","#E69F00"))+
theme_classic() +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank()) +
theme(legend.text = element_text(face = "italic"))+
stat_pvalue_manual(Post_hoc, hide.ns = TRUE) +
labs(
subtitle = get_test_label(kruskal.t, detailed = TRUE),
caption = get_pwc_label(Post_hoc))
Generalized Linear Model with quasi-Poisson distribution (q = 1)
M1 <- glm(q1 ~ SEASON*SPECIES+ELEVATION+SVL+SEQDEPTH,
family = quasipoisson(link = "log"),
data = Micro_div)
summary(M1)##
## Call:
## glm(formula = q1 ~ SEASON * SPECIES + ELEVATION + SVL + SEQDEPTH,
## family = quasipoisson(link = "log"), data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.262828 0.228362 18.667 < 2e-16
## SEASONRainy 0.405260 0.151684 2.672 0.009061
## SPECIESSceloporus bicanthalis 0.989264 0.459589 2.152 0.034227
## SPECIESSceloporus grammicus 0.028572 0.168483 0.170 0.865747
## SPECIESSceloporus spinosus 0.044405 0.217178 0.204 0.838486
## ELEVATION -0.387007 0.399665 -0.968 0.335661
## SVL 0.002237 0.003986 0.561 0.576205
## SEQDEPTH 0.235518 0.065659 3.587 0.000561
## SEASONRainy:SPECIESSceloporus bicanthalis -0.641160 0.190178 -3.371 0.001132
## SEASONRainy:SPECIESSceloporus grammicus -0.249934 0.196457 -1.272 0.206811
## SEASONRainy:SPECIESSceloporus spinosus -0.239895 0.202098 -1.187 0.238567
##
## (Intercept) ***
## SEASONRainy **
## SPECIESSceloporus bicanthalis *
## SPECIESSceloporus grammicus
## SPECIESSceloporus spinosus
## ELEVATION
## SVL
## SEQDEPTH ***
## SEASONRainy:SPECIESSceloporus bicanthalis **
## SEASONRainy:SPECIESSceloporus grammicus
## SEASONRainy:SPECIESSceloporus spinosus
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 11.08586)
##
## Null deviance: 1537.22 on 94 degrees of freedom
## Residual deviance: 999.94 on 84 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + ELEVATION + SVL + SEQDEPTH,
## family = quasipoisson(link = "log"), data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.462617 0.218130 20.458 < 2e-16 ***
## SEASONRainy 0.114355 0.078245 1.462 0.14748
## SPECIESSceloporus bicanthalis 0.652287 0.476593 1.369 0.17463
## SPECIESSceloporus grammicus -0.144594 0.116591 -1.240 0.21824
## SPECIESSceloporus spinosus -0.114180 0.186455 -0.612 0.54189
## ELEVATION -0.467503 0.427006 -1.095 0.27661
## SVL 0.001759 0.004152 0.424 0.67284
## SEQDEPTH 0.205901 0.068543 3.004 0.00348 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 12.24543)
##
## Null deviance: 1537.2 on 94 degrees of freedom
## Residual deviance: 1136.0 on 87 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + ELEVATION + SEQDEPTH, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.53327 0.14023 32.327 < 2e-16 ***
## SEASONRainy 0.11649 0.07772 1.499 0.13747
## SPECIESSceloporus bicanthalis 0.66616 0.47539 1.401 0.16464
## SPECIESSceloporus grammicus -0.11718 0.09653 -1.214 0.22803
## SPECIESSceloporus spinosus -0.04897 0.10416 -0.470 0.63943
## ELEVATION -0.48877 0.42406 -1.153 0.25220
## SEQDEPTH 0.20354 0.06796 2.995 0.00356 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 12.13221)
##
## Null deviance: 1537.2 on 94 degrees of freedom
## Residual deviance: 1138.2 on 88 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + SEQDEPTH, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.66137 0.08481 54.964 <2e-16 ***
## SEASONRainy 0.11993 0.07821 1.533 0.1287
## SPECIESSceloporus bicanthalis 0.12860 0.09341 1.377 0.1721
## SPECIESSceloporus grammicus -0.11600 0.09699 -1.196 0.2349
## SPECIESSceloporus spinosus -0.01756 0.10109 -0.174 0.8625
## SEQDEPTH 0.20435 0.06833 2.991 0.0036 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 12.2577)
##
## Null deviance: 1537.2 on 94 degrees of freedom
## Residual deviance: 1157.3 on 89 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 4.49221597 4.82475124
## SEASONRainy -0.03336187 0.27330751
## SPECIESSceloporus bicanthalis -0.05390201 0.31251946
## SPECIESSceloporus grammicus -0.30596309 0.07452011
## SPECIESSceloporus spinosus -0.21612307 0.18046386
## SEQDEPTH 0.06781350 0.33576942##
## Call:
## glm(formula = q1 ~ SEASON + SPECIES, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.60321 0.08660 53.155 < 2e-16 ***
## SEASONRainy 0.22149 0.07245 3.057 0.00294 **
## SPECIESSceloporus bicanthalis 0.15745 0.09734 1.618 0.10927
## SPECIESSceloporus grammicus -0.11983 0.10129 -1.183 0.23993
## SPECIESSceloporus spinosus -0.01096 0.10574 -0.104 0.91766
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 13.41224)
##
## Null deviance: 1537.2 on 94 degrees of freedom
## Residual deviance: 1260.8 on 90 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Best Model
M6 <- glm(q1 ~ SEASON*SPECIES,
family = quasipoisson(link = "log"),
data = Micro_div)
anova(M6, test="F")##
## Call:
## glm(formula = q1 ~ SEASON * SPECIES, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3930 0.1316 33.371 < 2e-16
## SEASONRainy 0.5318 0.1561 3.406 0.00100
## SPECIESSceloporus bicanthalis 0.5414 0.1656 3.270 0.00154
## SPECIESSceloporus grammicus 0.1005 0.1635 0.615 0.54032
## SPECIESSceloporus spinosus 0.1593 0.1715 0.929 0.35550
## SEASONRainy:SPECIESSceloporus bicanthalis -0.5796 0.2020 -2.869 0.00517
## SEASONRainy:SPECIESSceloporus grammicus -0.3298 0.2065 -1.597 0.11392
## SEASONRainy:SPECIESSceloporus spinosus -0.2364 0.2153 -1.098 0.27537
##
## (Intercept) ***
## SEASONRainy **
## SPECIESSceloporus bicanthalis **
## SPECIESSceloporus grammicus
## SPECIESSceloporus spinosus
## SEASONRainy:SPECIESSceloporus bicanthalis **
## SEASONRainy:SPECIESSceloporus grammicus
## SEASONRainy:SPECIESSceloporus spinosus
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 12.61408)
##
## Null deviance: 1537.2 on 94 degrees of freedom
## Residual deviance: 1151.1 on 87 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 4.1233396 4.64031820
## SEASONRainy 0.2317199 0.84503460
## SPECIESSceloporus bicanthalis 0.2207910 0.87128269
## SPECIESSceloporus grammicus -0.2155594 0.42677045
## SPECIESSceloporus spinosus -0.1742816 0.49967406
## SEASONRainy:SPECIESSceloporus bicanthalis -0.9789785 -0.18613994
## SEASONRainy:SPECIESSceloporus grammicus -0.7385988 0.07179439
## SEASONRainy:SPECIESSceloporus spinosus -0.6616867 0.18335432## (Intercept)
## 4.3929537
## SEASONRainy
## 0.5317728
## SPECIESSceloporus bicanthalis
## 0.5413804
## SPECIESSceloporus grammicus
## 0.1005116
## SPECIESSceloporus spinosus
## 0.1593282
## SEASONRainy:SPECIESSceloporus bicanthalis
## -0.5795744
## SEASONRainy:SPECIESSceloporus grammicus
## -0.3297727
## SEASONRainy:SPECIESSceloporus spinosus
## -0.2363617Residual plot of the best model
Plot model assumption
residuals_q1 <- plot_model(M6, type = "eff", terms = c("SEASON","SPECIES"),
colors = c("#999999", "#E69F00", "#56B4E9",
"#009E73", "#F0E442")) +
theme_classic() +
theme(legend.text = element_text(face = "italic")) +
theme(axis.text.x = element_text(size = 11))
residuals_q1 
Contrasts (dry vs rainy), CI95%, p values
Tukey Test
## SPECIES SEASON rate SE df asymp.LCL asymp.UCL
## Sceloporus aeneus Dry 80.9 10.65 Inf 62.5 105
## Sceloporus bicanthalis Dry 139.0 13.96 Inf 114.1 169
## Sceloporus grammicus Dry 89.4 8.67 Inf 74.0 108
## Sceloporus spinosus Dry 94.8 10.43 Inf 76.5 118
## Sceloporus aeneus Rainy 137.7 11.56 Inf 116.8 162
## Sceloporus bicanthalis Rainy 132.5 10.56 Inf 113.3 155
## Sceloporus grammicus Rainy 109.4 10.31 Inf 91.0 132
## Sceloporus spinosus Rainy 127.4 12.68 Inf 104.9 155
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaleM1.emm.s <- emmeans(M6, specs = ~ SPECIES*SEASON)
pairs(M1.emm.s, adjust = "tukey", infer=c(TRUE,TRUE))## contrast estimate SE df
## Sceloporus aeneus Dry - Sceloporus bicanthalis Dry -0.54138 0.166 Inf
## Sceloporus aeneus Dry - Sceloporus grammicus Dry -0.10051 0.164 Inf
## Sceloporus aeneus Dry - Sceloporus spinosus Dry -0.15933 0.172 Inf
## Sceloporus aeneus Dry - Sceloporus aeneus Rainy -0.53177 0.156 Inf
## Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy -0.49358 0.154 Inf
## Sceloporus aeneus Dry - Sceloporus grammicus Rainy -0.30251 0.162 Inf
## Sceloporus aeneus Dry - Sceloporus spinosus Rainy -0.45474 0.165 Inf
## Sceloporus bicanthalis Dry - Sceloporus grammicus Dry 0.44087 0.140 Inf
## Sceloporus bicanthalis Dry - Sceloporus spinosus Dry 0.38205 0.149 Inf
## Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy 0.00961 0.131 Inf
## Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy 0.04780 0.128 Inf
## Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy 0.23887 0.138 Inf
## Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy 0.08664 0.141 Inf
## Sceloporus grammicus Dry - Sceloporus spinosus Dry -0.05882 0.147 Inf
## Sceloporus grammicus Dry - Sceloporus aeneus Rainy -0.43126 0.128 Inf
## Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy -0.39307 0.125 Inf
## Sceloporus grammicus Dry - Sceloporus grammicus Rainy -0.20200 0.135 Inf
## Sceloporus grammicus Dry - Sceloporus spinosus Rainy -0.35423 0.139 Inf
## Sceloporus spinosus Dry - Sceloporus aeneus Rainy -0.37244 0.138 Inf
## Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy -0.33425 0.136 Inf
## Sceloporus spinosus Dry - Sceloporus grammicus Rainy -0.14318 0.145 Inf
## Sceloporus spinosus Dry - Sceloporus spinosus Rainy -0.29541 0.148 Inf
## Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy 0.03819 0.116 Inf
## Sceloporus aeneus Rainy - Sceloporus grammicus Rainy 0.22926 0.126 Inf
## Sceloporus aeneus Rainy - Sceloporus spinosus Rainy 0.07703 0.130 Inf
## Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy 0.19107 0.123 Inf
## Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy 0.03884 0.127 Inf
## Sceloporus grammicus Rainy - Sceloporus spinosus Rainy -0.15223 0.137 Inf
## asymp.LCL asymp.UCL z.ratio p.value
## -1.0432 -0.0396 -3.270 0.0239
## -0.5961 0.3950 -0.615 0.9987
## -0.6792 0.3605 -0.929 0.9833
## -1.0050 -0.0585 -3.406 0.0152
## -0.9599 -0.0272 -3.208 0.0292
## -0.7931 0.1880 -1.869 0.5723
## -0.9548 0.0454 -2.756 0.1064
## 0.0178 0.8640 3.158 0.0341
## -0.0693 0.8334 2.566 0.1686
## -0.3871 0.4063 0.073 1.0000
## -0.3407 0.4363 0.373 1.0000
## -0.1784 0.6561 1.735 0.6640
## -0.3418 0.5151 0.613 0.9987
## -0.5032 0.3855 -0.401 0.9999
## -0.8200 -0.0425 -3.362 0.0176
## -0.7734 -0.0127 -3.132 0.0369
## -0.6117 0.2077 -1.495 0.8107
## -0.7753 0.0668 -2.550 0.1748
## -0.7917 0.0469 -2.692 0.1248
## -0.7458 0.0773 -2.462 0.2121
## -0.5819 0.2956 -0.989 0.9761
## -0.7448 0.1540 -1.992 0.4871
## -0.3126 0.3890 0.330 1.0000
## -0.1531 0.6116 1.817 0.6082
## -0.3175 0.4716 0.592 0.9990
## -0.1828 0.5649 1.549 0.7805
## -0.3475 0.4251 0.305 1.0000
## -0.5674 0.2629 -1.111 0.9546
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 8 estimates
## P value adjustment: tukey method for comparing a family of 8 estimates## Sceloporus aeneus Dry Sceloporus bicanthalis Dry
## Sceloporus aeneus Dry [4.39] 0.0239
## Sceloporus bicanthalis Dry -0.54138 [4.93]
## Sceloporus grammicus Dry -0.10051 0.44087
## Sceloporus spinosus Dry -0.15933 0.38205
## Sceloporus aeneus Rainy -0.53177 0.00961
## Sceloporus bicanthalis Rainy -0.49358 0.04780
## Sceloporus grammicus Rainy -0.30251 0.23887
## Sceloporus spinosus Rainy -0.45474 0.08664
## Sceloporus grammicus Dry Sceloporus spinosus Dry
## Sceloporus aeneus Dry 0.9987 0.9833
## Sceloporus bicanthalis Dry 0.0341 0.1686
## Sceloporus grammicus Dry [4.49] 0.9999
## Sceloporus spinosus Dry -0.05882 [4.55]
## Sceloporus aeneus Rainy -0.43126 -0.37244
## Sceloporus bicanthalis Rainy -0.39307 -0.33425
## Sceloporus grammicus Rainy -0.20200 -0.14318
## Sceloporus spinosus Rainy -0.35423 -0.29541
## Sceloporus aeneus Rainy
## Sceloporus aeneus Dry 0.0152
## Sceloporus bicanthalis Dry 1.0000
## Sceloporus grammicus Dry 0.0176
## Sceloporus spinosus Dry 0.1248
## Sceloporus aeneus Rainy [4.92]
## Sceloporus bicanthalis Rainy 0.03819
## Sceloporus grammicus Rainy 0.22926
## Sceloporus spinosus Rainy 0.07703
## Sceloporus bicanthalis Rainy
## Sceloporus aeneus Dry 0.0292
## Sceloporus bicanthalis Dry 1.0000
## Sceloporus grammicus Dry 0.0369
## Sceloporus spinosus Dry 0.2121
## Sceloporus aeneus Rainy 1.0000
## Sceloporus bicanthalis Rainy [4.89]
## Sceloporus grammicus Rainy 0.19107
## Sceloporus spinosus Rainy 0.03884
## Sceloporus grammicus Rainy
## Sceloporus aeneus Dry 0.5723
## Sceloporus bicanthalis Dry 0.6640
## Sceloporus grammicus Dry 0.8107
## Sceloporus spinosus Dry 0.9761
## Sceloporus aeneus Rainy 0.6082
## Sceloporus bicanthalis Rainy 0.7805
## Sceloporus grammicus Rainy [4.70]
## Sceloporus spinosus Rainy -0.15223
## Sceloporus spinosus Rainy
## Sceloporus aeneus Dry 0.1064
## Sceloporus bicanthalis Dry 0.9987
## Sceloporus grammicus Dry 0.1748
## Sceloporus spinosus Dry 0.4871
## Sceloporus aeneus Rainy 0.9990
## Sceloporus bicanthalis Rainy 1.0000
## Sceloporus grammicus Rainy 0.9546
## Sceloporus spinosus Rainy [4.85]
##
## Row and column labels: SPECIES:SEASON
## Upper triangle: P values adjust = "tukey"
## Diagonal: [Estimates] (emmean)
## Lower triangle: Comparisons (estimate) earlier vs. later
Evaluate whether sex variable influence on gut microbiota diversity (q=2)
# Paired test (Wilcoxon test, q2)
# Sceloporus aeneus
aeneus_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" &
Sex == "Female", q2, drop = TRUE)
aeneus_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" &
Sex == "Male", q2, drop = TRUE)
aeneus_FemMale_GutMic <- wilcox.test(x= aeneus_fem_GutMic, y= aeneus_male_GutMic)
aeneus_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: aeneus_fem_GutMic and aeneus_male_GutMic
## W = 39, p-value = 0.1713
## alternative hypothesis: true location shift is not equal to 0# Sceloporus bicanthalis
bica_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" &
Sex == "Female", q2, drop = TRUE)
bica_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" &
Sex == "Male", q2, drop = TRUE)
bica_FemMale_GutMic <- wilcox.test(x= bica_fem_GutMic, y= bica_male_GutMic)
bica_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: bica_fem_GutMic and bica_male_GutMic
## W = 69, p-value = 0.9095
## alternative hypothesis: true location shift is not equal to 0# Sceloporus grammicus
gram_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" &
Sex == "Female", q2, drop = TRUE)
gram_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" &
Sex == "Male", q2, drop = TRUE)
gram_FemMale_GutMic <- wilcox.test(x= gram_fem_GutMic, y= gram_male_GutMic)
gram_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: gram_fem_GutMic and gram_male_GutMic
## W = 107, p-value = 0.6313
## alternative hypothesis: true location shift is not equal to 0# Sceloporus spinosus
spi_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" &
Sex == "Female", q2, drop = TRUE)
spi_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" &
Sex == "Male", q2, drop = TRUE)
spi_FemMale_GutMic <- wilcox.test(x= spi_fem_GutMic, y= spi_male_GutMic)
spi_FemMale_GutMic##
## Wilcoxon rank sum exact test
##
## data: spi_fem_GutMic and spi_male_GutMic
## W = 40, p-value = 0.4137
## alternative hypothesis: true location shift is not equal to 0Generalized Linear Model with quasi-Poisson distribution (q = 2)
MODEL1 <- glm(q2 ~ SEASON*SPECIES+ELEVATION+SVL+SEQDEPTH,
family = quasipoisson(link = "log"),
data = Micro_div)
summary(MODEL1)##
## Call:
## glm(formula = q2 ~ SEASON * SPECIES + ELEVATION + SVL + SEQDEPTH,
## family = quasipoisson(link = "log"), data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3.489e+00 3.263e-01 10.694
## SEASONRainy 3.432e-01 2.143e-01 1.601
## SPECIESSceloporus bicanthalis 1.267e+00 6.553e-01 1.933
## SPECIESSceloporus grammicus 6.931e-02 2.358e-01 0.294
## SPECIESSceloporus spinosus 1.270e-01 3.062e-01 0.415
## ELEVATION -4.765e-01 5.739e-01 -0.830
## SVL -2.387e-05 5.754e-03 -0.004
## SEQDEPTH 1.204e-01 1.009e-01 1.193
## SEASONRainy:SPECIESSceloporus bicanthalis -8.202e-01 2.652e-01 -3.093
## SEASONRainy:SPECIESSceloporus grammicus -1.219e-01 2.754e-01 -0.442
## SEASONRainy:SPECIESSceloporus spinosus -3.317e-01 2.909e-01 -1.140
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## SEASONRainy 0.11305
## SPECIESSceloporus bicanthalis 0.05657 .
## SPECIESSceloporus grammicus 0.76956
## SPECIESSceloporus spinosus 0.67944
## ELEVATION 0.40867
## SVL 0.99670
## SEQDEPTH 0.23636
## SEASONRainy:SPECIESSceloporus bicanthalis 0.00269 **
## SEASONRainy:SPECIESSceloporus grammicus 0.65931
## SEASONRainy:SPECIESSceloporus spinosus 0.25746
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 9.481509)
##
## Null deviance: 1039.03 on 94 degrees of freedom
## Residual deviance: 818.79 on 84 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION + SVL + SEQDEPTH,
## family = quasipoisson(link = "log"), data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.708501 0.314372 11.797 <2e-16 ***
## SEASONRainy 0.025511 0.109991 0.232 0.817
## SPECIESSceloporus bicanthalis 0.933649 0.676327 1.380 0.171
## SPECIESSceloporus grammicus -0.025499 0.165697 -0.154 0.878
## SPECIESSceloporus spinosus -0.051587 0.267991 -0.192 0.848
## ELEVATION -0.626981 0.606572 -1.034 0.304
## SVL -0.001225 0.005989 -0.205 0.838
## SEQDEPTH 0.041928 0.106778 0.393 0.696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.47637)
##
## Null deviance: 1039.03 on 94 degrees of freedom
## Residual deviance: 934.33 on 87 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION + SEQDEPTH, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.65897 0.19913 18.375 <2e-16 ***
## SEASONRainy 0.02422 0.10921 0.222 0.825
## SPECIESSceloporus bicanthalis 0.92494 0.66933 1.382 0.171
## SPECIESSceloporus grammicus -0.04438 0.13685 -0.324 0.746
## SPECIESSceloporus spinosus -0.09673 0.15179 -0.637 0.526
## ELEVATION -0.61281 0.59747 -1.026 0.308
## SEQDEPTH 0.04319 0.10604 0.407 0.685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.36163)
##
## Null deviance: 1039.03 on 94 degrees of freedom
## Residual deviance: 934.77 on 88 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.64561 0.19533 18.664 <2e-16 ***
## SEASONRainy 0.04488 0.09589 0.468 0.641
## SPECIESSceloporus bicanthalis 0.93231 0.66505 1.402 0.164
## SPECIESSceloporus grammicus -0.04395 0.13615 -0.323 0.748
## SPECIESSceloporus spinosus -0.09462 0.15095 -0.627 0.532
## ELEVATION -0.61383 0.59391 -1.034 0.304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.26101)
##
## Null deviance: 1039.03 on 94 degrees of freedom
## Residual deviance: 936.46 on 89 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SEASON + SPECIES, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.80662 0.11637 32.712 <2e-16 ***
## SEASONRainy 0.04940 0.09642 0.512 0.6097
## SPECIESSceloporus bicanthalis 0.25823 0.13175 1.960 0.0531 .
## SPECIESSceloporus grammicus -0.04288 0.13674 -0.314 0.7546
## SPECIESSceloporus spinosus -0.05584 0.14700 -0.380 0.7049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.3626)
##
## Null deviance: 1039.03 on 94 degrees of freedom
## Residual deviance: 950.07 on 90 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Best Model
MODEL6 <- glm(q2 ~ SEASON*SPECIES,
family = quasipoisson(link = "log"),
data = Micro_div)
summary(MODEL6)##
## Call:
## glm(formula = q2 ~ SEASON * SPECIES, family = quasipoisson(link = "log"),
## data = Micro_div)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.57869 0.17122 20.902 < 2e-16
## SEASONRainy 0.40334 0.20706 1.948 0.054645
## SPECIESSceloporus bicanthalis 0.74146 0.20804 3.564 0.000596
## SPECIESSceloporus grammicus 0.08107 0.21339 0.380 0.704917
## SPECIESSceloporus spinosus 0.15897 0.22310 0.713 0.478026
## SEASONRainy:SPECIESSceloporus bicanthalis -0.79089 0.26304 -3.007 0.003451
## SEASONRainy:SPECIESSceloporus grammicus -0.14701 0.27125 -0.542 0.589223
## SEASONRainy:SPECIESSceloporus spinosus -0.32729 0.29016 -1.128 0.262430
##
## (Intercept) ***
## SEASONRainy .
## SPECIESSceloporus bicanthalis ***
## SPECIESSceloporus grammicus
## SPECIESSceloporus spinosus
## SEASONRainy:SPECIESSceloporus bicanthalis **
## SEASONRainy:SPECIESSceloporus grammicus
## SEASONRainy:SPECIESSceloporus spinosus
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 9.452614)
##
## Null deviance: 1039.03 on 94 degrees of freedom
## Residual deviance: 840.53 on 87 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4## (Intercept)
## 3.57869413
## SEASONRainy
## 0.40334451
## SPECIESSceloporus bicanthalis
## 0.74145835
## SPECIESSceloporus grammicus
## 0.08107416
## SPECIESSceloporus spinosus
## 0.15897437
## SEASONRainy:SPECIESSceloporus bicanthalis
## -0.79088870
## SEASONRainy:SPECIESSceloporus grammicus
## -0.14700948
## SEASONRainy:SPECIESSceloporus spinosus
## -0.32729337## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 3.223230360 3.8964819
## SEASONRainy 0.006183591 0.8207875
## SPECIESSceloporus bicanthalis 0.342095601 1.1605800
## SPECIESSceloporus grammicus -0.330307974 0.5093207
## SPECIESSceloporus spinosus -0.274369614 0.6037278
## SEASONRainy:SPECIESSceloporus bicanthalis -1.313319633 -0.2802605
## SEASONRainy:SPECIESSceloporus grammicus -0.685102597 0.3801804
## SEASONRainy:SPECIESSceloporus spinosus -0.902353379 0.2372010Residual plot of the best model
Plot model assumption
residuals_q2 <- plot_model(MODEL6, type = "eff", terms = c("SEASON","SPECIES"),
colors = c("#999999", "#E69F00", "#56B4E9",
"#009E73", "#F0E442")) +
theme_classic() +
theme(legend.text = element_text(face = "italic")) +
theme(axis.text.x = element_text(size = 11))
residuals_q2 
Contrasts (dry vs rainy), CI95%, p values
## SPECIES SEASON rate SE df asymp.LCL asymp.UCL
## Sceloporus aeneus Dry 35.8 6.13 Inf 25.6 50.1
## Sceloporus bicanthalis Dry 75.2 8.89 Inf 59.7 94.8
## Sceloporus grammicus Dry 38.9 4.95 Inf 30.3 49.9
## Sceloporus spinosus Dry 42.0 6.01 Inf 31.7 55.6
## Sceloporus aeneus Rainy 53.6 6.24 Inf 42.7 67.4
## Sceloporus bicanthalis Rainy 51.0 5.67 Inf 41.1 63.5
## Sceloporus grammicus Rainy 50.2 6.04 Inf 39.7 63.6
## Sceloporus spinosus Rainy 45.3 6.55 Inf 34.1 60.1
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaleM1.emm.s_q2 <- emmeans(MODEL6, specs = ~ SPECIES*SEASON)
pairs(M1.emm.s_q2, adjust = "tukey", infer=c(TRUE,TRUE))## contrast estimate SE df
## Sceloporus aeneus Dry - Sceloporus bicanthalis Dry -0.7415 0.208 Inf
## Sceloporus aeneus Dry - Sceloporus grammicus Dry -0.0811 0.213 Inf
## Sceloporus aeneus Dry - Sceloporus spinosus Dry -0.1590 0.223 Inf
## Sceloporus aeneus Dry - Sceloporus aeneus Rainy -0.4033 0.207 Inf
## Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy -0.3539 0.204 Inf
## Sceloporus aeneus Dry - Sceloporus grammicus Rainy -0.3374 0.209 Inf
## Sceloporus aeneus Dry - Sceloporus spinosus Rainy -0.2350 0.224 Inf
## Sceloporus bicanthalis Dry - Sceloporus grammicus Dry 0.6604 0.174 Inf
## Sceloporus bicanthalis Dry - Sceloporus spinosus Dry 0.5825 0.186 Inf
## Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy 0.3381 0.166 Inf
## Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy 0.3875 0.162 Inf
## Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy 0.4040 0.169 Inf
## Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy 0.5064 0.187 Inf
## Sceloporus grammicus Dry - Sceloporus spinosus Dry -0.0779 0.192 Inf
## Sceloporus grammicus Dry - Sceloporus aeneus Rainy -0.3223 0.173 Inf
## Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy -0.2728 0.169 Inf
## Sceloporus grammicus Dry - Sceloporus grammicus Rainy -0.2563 0.175 Inf
## Sceloporus grammicus Dry - Sceloporus spinosus Rainy -0.1540 0.193 Inf
## Sceloporus spinosus Dry - Sceloporus aeneus Rainy -0.2444 0.184 Inf
## Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy -0.1949 0.181 Inf
## Sceloporus spinosus Dry - Sceloporus grammicus Rainy -0.1784 0.187 Inf
## Sceloporus spinosus Dry - Sceloporus spinosus Rainy -0.0761 0.203 Inf
## Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy 0.0494 0.161 Inf
## Sceloporus aeneus Rainy - Sceloporus grammicus Rainy 0.0659 0.167 Inf
## Sceloporus aeneus Rainy - Sceloporus spinosus Rainy 0.1683 0.186 Inf
## Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy 0.0165 0.164 Inf
## Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy 0.1189 0.182 Inf
## Sceloporus grammicus Rainy - Sceloporus spinosus Rainy 0.1024 0.188 Inf
## asymp.LCL asymp.UCL z.ratio p.value
## -1.3720 -0.111 -3.564 0.0087
## -0.7278 0.566 -0.380 0.9999
## -0.8352 0.517 -0.713 0.9967
## -1.0309 0.224 -1.948 0.5176
## -0.9726 0.265 -1.734 0.6649
## -0.9717 0.297 -1.612 0.7431
## -0.9139 0.444 -1.049 0.9668
## 0.1338 1.187 3.801 0.0036
## 0.0201 1.145 3.139 0.0361
## -0.1647 0.841 2.038 0.4560
## -0.1041 0.879 2.389 0.2465
## -0.1072 0.915 2.396 0.2433
## -0.0592 1.072 2.714 0.1183
## -0.6584 0.503 -0.407 0.9999
## -0.8453 0.201 -1.868 0.5734
## -0.7851 0.239 -1.614 0.7418
## -0.7874 0.275 -1.463 0.8272
## -0.7376 0.430 -0.800 0.9932
## -0.8034 0.315 -1.325 0.8897
## -0.7439 0.354 -1.076 0.9618
## -0.7450 0.388 -0.955 0.9805
## -0.6921 0.540 -0.374 1.0000
## -0.4384 0.537 0.307 1.0000
## -0.4416 0.573 0.394 0.9999
## -0.3940 0.731 0.907 0.9855
## -0.4799 0.513 0.101 1.0000
## -0.4334 0.671 0.652 0.9981
## -0.4674 0.672 0.545 0.9994
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 8 estimates
## P value adjustment: tukey method for comparing a family of 8 estimates## contrast estimate SE df
## Sceloporus aeneus Dry - Sceloporus bicanthalis Dry -0.7415 0.208 Inf
## Sceloporus aeneus Dry - Sceloporus grammicus Dry -0.0811 0.213 Inf
## Sceloporus aeneus Dry - Sceloporus spinosus Dry -0.1590 0.223 Inf
## Sceloporus aeneus Dry - Sceloporus aeneus Rainy -0.4033 0.207 Inf
## Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy -0.3539 0.204 Inf
## Sceloporus aeneus Dry - Sceloporus grammicus Rainy -0.3374 0.209 Inf
## Sceloporus aeneus Dry - Sceloporus spinosus Rainy -0.2350 0.224 Inf
## Sceloporus bicanthalis Dry - Sceloporus grammicus Dry 0.6604 0.174 Inf
## Sceloporus bicanthalis Dry - Sceloporus spinosus Dry 0.5825 0.186 Inf
## Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy 0.3381 0.166 Inf
## Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy 0.3875 0.162 Inf
## Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy 0.4040 0.169 Inf
## Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy 0.5064 0.187 Inf
## Sceloporus grammicus Dry - Sceloporus spinosus Dry -0.0779 0.192 Inf
## Sceloporus grammicus Dry - Sceloporus aeneus Rainy -0.3223 0.173 Inf
## Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy -0.2728 0.169 Inf
## Sceloporus grammicus Dry - Sceloporus grammicus Rainy -0.2563 0.175 Inf
## Sceloporus grammicus Dry - Sceloporus spinosus Rainy -0.1540 0.193 Inf
## Sceloporus spinosus Dry - Sceloporus aeneus Rainy -0.2444 0.184 Inf
## Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy -0.1949 0.181 Inf
## Sceloporus spinosus Dry - Sceloporus grammicus Rainy -0.1784 0.187 Inf
## Sceloporus spinosus Dry - Sceloporus spinosus Rainy -0.0761 0.203 Inf
## Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy 0.0494 0.161 Inf
## Sceloporus aeneus Rainy - Sceloporus grammicus Rainy 0.0659 0.167 Inf
## Sceloporus aeneus Rainy - Sceloporus spinosus Rainy 0.1683 0.186 Inf
## Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy 0.0165 0.164 Inf
## Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy 0.1189 0.182 Inf
## Sceloporus grammicus Rainy - Sceloporus spinosus Rainy 0.1024 0.188 Inf
## z.ratio p.value
## -3.564 0.0087
## -0.380 0.9999
## -0.713 0.9967
## -1.948 0.5176
## -1.734 0.6649
## -1.612 0.7431
## -1.049 0.9668
## 3.801 0.0036
## 3.139 0.0361
## 2.038 0.4560
## 2.389 0.2465
## 2.396 0.2433
## 2.714 0.1183
## -0.407 0.9999
## -1.868 0.5734
## -1.614 0.7418
## -1.463 0.8272
## -0.800 0.9932
## -1.325 0.8897
## -1.076 0.9618
## -0.955 0.9805
## -0.374 1.0000
## 0.307 1.0000
## 0.394 0.9999
## 0.907 0.9855
## 0.101 1.0000
## 0.652 0.9981
## 0.545 0.9994
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 8 estimates## Sceloporus aeneus Dry Sceloporus bicanthalis Dry
## Sceloporus aeneus Dry [3.58] 0.0087
## Sceloporus bicanthalis Dry -0.7415 [4.32]
## Sceloporus grammicus Dry -0.0811 0.6604
## Sceloporus spinosus Dry -0.1590 0.5825
## Sceloporus aeneus Rainy -0.4033 0.3381
## Sceloporus bicanthalis Rainy -0.3539 0.3875
## Sceloporus grammicus Rainy -0.3374 0.4040
## Sceloporus spinosus Rainy -0.2350 0.5064
## Sceloporus grammicus Dry Sceloporus spinosus Dry
## Sceloporus aeneus Dry 0.9999 0.9967
## Sceloporus bicanthalis Dry 0.0036 0.0361
## Sceloporus grammicus Dry [3.66] 0.9999
## Sceloporus spinosus Dry -0.0779 [3.74]
## Sceloporus aeneus Rainy -0.3223 -0.2444
## Sceloporus bicanthalis Rainy -0.2728 -0.1949
## Sceloporus grammicus Rainy -0.2563 -0.1784
## Sceloporus spinosus Rainy -0.1540 -0.0761
## Sceloporus aeneus Rainy
## Sceloporus aeneus Dry 0.5176
## Sceloporus bicanthalis Dry 0.4560
## Sceloporus grammicus Dry 0.5734
## Sceloporus spinosus Dry 0.8897
## Sceloporus aeneus Rainy [3.98]
## Sceloporus bicanthalis Rainy 0.0494
## Sceloporus grammicus Rainy 0.0659
## Sceloporus spinosus Rainy 0.1683
## Sceloporus bicanthalis Rainy
## Sceloporus aeneus Dry 0.6649
## Sceloporus bicanthalis Dry 0.2465
## Sceloporus grammicus Dry 0.7418
## Sceloporus spinosus Dry 0.9618
## Sceloporus aeneus Rainy 1.0000
## Sceloporus bicanthalis Rainy [3.93]
## Sceloporus grammicus Rainy 0.0165
## Sceloporus spinosus Rainy 0.1189
## Sceloporus grammicus Rainy
## Sceloporus aeneus Dry 0.7431
## Sceloporus bicanthalis Dry 0.2433
## Sceloporus grammicus Dry 0.8272
## Sceloporus spinosus Dry 0.9805
## Sceloporus aeneus Rainy 0.9999
## Sceloporus bicanthalis Rainy 1.0000
## Sceloporus grammicus Rainy [3.92]
## Sceloporus spinosus Rainy 0.1024
## Sceloporus spinosus Rainy
## Sceloporus aeneus Dry 0.9668
## Sceloporus bicanthalis Dry 0.1183
## Sceloporus grammicus Dry 0.9932
## Sceloporus spinosus Dry 1.0000
## Sceloporus aeneus Rainy 0.9855
## Sceloporus bicanthalis Rainy 0.9981
## Sceloporus grammicus Rainy 0.9994
## Sceloporus spinosus Rainy [3.81]
##
## Row and column labels: SPECIES:SEASON
## Upper triangle: P values adjust = "tukey"
## Diagonal: [Estimates] (emmean)
## Lower triangle: Comparisons (estimate) earlier vs. later
Relation between alpha taxonomic microbiota and diet
Micro_div$Ind <- substr(Micro_div$SampleID, 8, nchar(Micro_div$SampleID) - 0)
diet_phylo_div <- read.csv("../data/diet_Hill_numbers_q012.csv", header = TRUE)
# Generar una nueva columna con la información de otra columna, excluyendo los últimos 4 caracteres
diet_phylo_div$Ind <- substr(diet_phylo_div$SampleID, 8, nchar(diet_phylo_div$SampleID) - 4)
#names(diet_phylo_div)[names(diet_phylo_div) %in% c("q0", "q1", "q2")] <- c("d.p.0", "d.p.1", #"d.p.2")
divs <- Micro_div %>%
inner_join(diet_phylo_div, by = c("Ind"="Ind"))
#write.table(divs, file="../data/full_divs.txt", sep = "\t")
##Plot
taxq1_tq0 <- ggscatter(divs, x = "tq0", y = "q1", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), xlab= "Taxonomic richness of diet (number of taxonomic groups)", ylab="Gut bacterial diversity (effective number of common ASVs)",
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")
taxq2_tq0 <- ggscatter(divs, x = "tq0", y = "q2", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"),
xlab= "Taxonomic richness of diet (number of taxonomic groups)", ylab="Gut bacterial diversity (effective number of dominant ASVs)",
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")
phyq1_pq0 <- ggscatter(divs, x = "pq0", y = "q1", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73") ,
xlab= " Phylogenetic richness of diet (number of phylogenetic lineages)", ylab="Gut bacterial diversity (effective number of common ASVs)",
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")
phyq2_pq0 <- ggscatter(divs, x = "pq0", y = "q2", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"),
xlab= "Phylogenetic richness of diet (number of phylogenetic lineages)",
ylab="Gut bacterial diversity (effective number of dominant ASVs)",
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")
taxq1_tq0.sp <- ggscatter(divs, x = "tq0", y = "q1", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"),
xlab= "Taxonomic richness of diet (number of taxonomic groups)",
ylab="Gut bacterial diversity (effective number of common ASVs)",
facet.by = "Species_Scel",
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson",
label.x = 3, label.sep = "\n")
)+
theme(legend.position = "none",
strip.text = element_text(face = "italic"))
taxq2_tq0.sp <- ggscatter(divs, x = "tq0", y = "q2", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"),
xlab= "Taxonomic richness of diet (number of taxonomic groups)", ylab="Gut bacterial diversity (effective number of dominant ASVs)",
facet.by = "Species_Scel" ,
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.position = "none",
strip.text = element_text(face = "italic"))
phyq1_pq0.sp <- ggscatter(divs, x = "pq0", y = "q1", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73") ,
xlab= " Phylogenetic richness of diet (number of phylogenetic lineages)", ylab="Gut bacterial diversity (effective number of common ASVs)",
facet.by = "Species_Scel",
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.position = "none",
strip.text = element_text(face = "italic"))
phyq2_pq0.sp <- ggscatter(divs, x = "pq0", y = "q2", color = "Species_Scel",
palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"),
xlab= "Phylogenetic richness of diet (number of phylogenetic lineages)", ylab="Gut bacterial diversity (effective number of dominant ASVs)",
facet.by = "Species_Scel" ,
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.position = "none",
strip.text = element_text(face = "italic"))
taxq1_tq0
diet.vs.micro1<- plot_grid(taxq1_tq0, taxq1_tq0.sp,
labels = "AUTO", ncol = 1)
#ggsave("../figures/dietvsmicro1.jpeg", width=8, height=11, dpi=300)
diet.vs.micro2<- plot_grid(taxq2_tq0, taxq2_tq0.sp,
labels = c("E","F"), ncol = 1)
#ggsave("../figures/dietvsmicro2.jpeg", width=8, height=11, dpi=300)
diet.vs.micro3<- plot_grid(phyq1_pq0, phyq1_pq0.sp,
labels = c("C", "D"), ncol = 1)
#ggsave("../figures/dietvsmicro3.jpeg", width=8, height=11, dpi=300)
diet.vs.micro4<- plot_grid(phyq2_pq0, phyq2_pq0.sp,
labels = c("G", "H"), ncol = 1)
#ggsave("../figures/dietvsmicro4.jpeg", width=8, height=11, dpi=300)Effect of diet on gut microbiota
loading files
Check distribution of the dependent variable (diversity q=0, Hill numbers)
##
## Shapiro-Wilk normality test
##
## data: Microb_data$q0
## W = 0.92649, p-value = 0.0009384##
## Shapiro-Wilk normality test
##
## data: Microb_data$q1
## W = 0.97555, p-value = 0.2334
##
## Bartlett test of homogeneity of variances
##
## data: q0 by Species_Scel
## Bartlett's K-squared = 9.0931, df = 3, p-value = 0.02808## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.Plot variances among lizard species (q0)
boxplot(q0 ~ Species_Scel, data = Microb_data,
main = "Differences in variance among lizard species",
xlab = "Species",
ylab = " Hill numbers (q0)",
col = "#4682B433",
border = "black")
###Preparing the data
#Declare as factors
SPECIES <- as.factor(Microb_data$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(Microb_data$Season) # Dry and Rainy
SEX <- as.factor(Microb_data$Sex) # Male and Female
#Standardize continous independent variables
SVL <- rescale(Microb_data$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(Microb_data$Elevation, binary.inputs = "center") # Taken as m a.s.l.
SEQDEPTH <- rescale(Microb_data$SeqDepth, binary.inputs = "center")Linear model with Poisson distribution (taxonomic q = 0)
With Microbial Diversity (q1)
## UPDATE FUNCTION (RUN ANOVA)
# 68 samples (diet / microbiota)
m1 <- glm(q1 ~ Diet*SPECIES+SEASON,
family = quasipoisson(link = "log"),
data = Microb_data)
summary(m1)##
## Call:
## glm(formula = q1 ~ Diet * SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6562151 0.1865713 24.957 < 2e-16 ***
## Diet -0.0004683 0.0196703 -0.024 0.98109
## SPECIESSceloporus bicanthalis 0.0595147 0.2216353 0.269 0.78930
## SPECIESSceloporus grammicus -0.1948473 0.2453897 -0.794 0.43059
## SPECIESSceloporus spinosus 0.1751032 0.2349154 0.745 0.45921
## SEASONRainy 0.2689010 0.0857259 3.137 0.00274 **
## Diet:SPECIESSceloporus bicanthalis 0.0018106 0.0258806 0.070 0.94448
## Diet:SPECIESSceloporus grammicus -0.0084125 0.0312962 -0.269 0.78909
## Diet:SPECIESSceloporus spinosus -0.0373851 0.0318303 -1.175 0.24525
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.57198)
##
## Null deviance: 932.15 on 63 degrees of freedom
## Residual deviance: 593.29 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = q1 ~ Diet + SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.704855 0.130358 36.092 < 2e-16 ***
## Diet -0.007877 0.010173 -0.774 0.44187
## SPECIESSceloporus bicanthalis 0.061506 0.106646 0.577 0.56636
## SPECIESSceloporus grammicus -0.256866 0.111237 -2.309 0.02452 *
## SPECIESSceloporus spinosus -0.047488 0.116892 -0.406 0.68605
## SEASONRainy 0.283586 0.080238 3.534 0.00081 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.42645)
##
## Null deviance: 932.15 on 63 degrees of freedom
## Residual deviance: 614.21 on 58 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.63657 0.09649 48.053 < 2e-16 ***
## SPECIESSceloporus bicanthalis 0.07079 0.10567 0.670 0.505539
## SPECIESSceloporus grammicus -0.24165 0.10918 -2.213 0.030761 *
## SPECIESSceloporus spinosus -0.02901 0.11425 -0.254 0.800452
## SEASONRainy 0.29232 0.07918 3.692 0.000488 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.37193)
##
## Null deviance: 932.15 on 63 degrees of freedom
## Residual deviance: 620.50 on 59 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.55181 0.06272 72.574 < 2e-16 ***
## SEASONRainy 0.33820 0.08005 4.225 7.99e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 11.18809)
##
## Null deviance: 932.15 on 63 degrees of freedom
## Residual deviance: 728.20 on 62 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SPECIES, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.82632 0.08770 55.033 <2e-16 ***
## SPECIESSceloporus bicanthalis 0.09084 0.11575 0.785 0.4356
## SPECIESSceloporus grammicus -0.29554 0.11870 -2.490 0.0156 *
## SPECIESSceloporus spinosus -0.08290 0.12430 -0.667 0.5074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 12.47322)
##
## Null deviance: 932.15 on 63 degrees of freedom
## Residual deviance: 764.28 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Plot model assumption
plot_model(m3, type = "eff", terms = "SPECIES",
colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
theme(legend.text = element_text(face = "italic"))
plot_model(m3, type = "eff", terms = c("SPECIES", "SEASON"),
colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
theme_classic() +
theme(legend.text = element_text(face = "italic")) +
theme(axis.text.x = element_text(size = 11))
With Microbial Diversity (q2)
model1 <- glm(q2 ~ Diet*SPECIES+SEASON,
family = quasipoisson(link = "log"),
data = Microb_data)
summary(model1)##
## Call:
## glm(formula = q2 ~ Diet * SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.62441 0.25085 14.449 <2e-16 ***
## Diet 0.01764 0.02568 0.687 0.4950
## SPECIESSceloporus bicanthalis 0.34411 0.29428 1.169 0.2473
## SPECIESSceloporus grammicus 0.08006 0.32076 0.250 0.8038
## SPECIESSceloporus spinosus 0.45688 0.31449 1.453 0.1520
## SEASONRainy 0.11588 0.10904 1.063 0.2925
## Diet:SPECIESSceloporus bicanthalis -0.01057 0.03311 -0.319 0.7507
## Diet:SPECIESSceloporus grammicus -0.03106 0.04019 -0.773 0.4429
## Diet:SPECIESSceloporus spinosus -0.07581 0.04204 -1.803 0.0768 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.255196)
##
## Null deviance: 528.4 on 63 degrees of freedom
## Residual deviance: 401.7 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = q2 ~ Diet + SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.79639 0.17233 22.030 <2e-16 ***
## Diet -0.00507 0.01324 -0.383 0.7030
## SPECIESSceloporus bicanthalis 0.24216 0.14164 1.710 0.0927 .
## SPECIESSceloporus grammicus -0.15074 0.14887 -1.013 0.3155
## SPECIESSceloporus spinosus -0.02013 0.15835 -0.127 0.8993
## SEASONRainy 0.13451 0.10404 1.293 0.2012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.364622)
##
## Null deviance: 528.40 on 63 degrees of freedom
## Residual deviance: 430.43 on 58 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.752299 0.127845 29.350 <2e-16 ***
## SPECIESSceloporus bicanthalis 0.247615 0.139789 1.771 0.0817 .
## SPECIESSceloporus grammicus -0.140309 0.145194 -0.966 0.3378
## SPECIESSceloporus spinosus -0.008757 0.154400 -0.057 0.9550
## SEASONRainy 0.140365 0.102034 1.376 0.1741
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.248433)
##
## Null deviance: 528.40 on 63 degrees of freedom
## Residual deviance: 431.52 on 59 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SPECIES, family = quasipoisson(link = "log"),
## data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.84098 0.10960 35.046 <2e-16 ***
## SPECIESSceloporus bicanthalis 0.25751 0.13984 1.841 0.0705 .
## SPECIESSceloporus grammicus -0.16641 0.14420 -1.154 0.2531
## SPECIESSceloporus spinosus -0.03486 0.15350 -0.227 0.8211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.272215)
##
## Null deviance: 528.40 on 63 degrees of freedom
## Residual deviance: 445.32 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ 1, family = quasipoisson(link = "log"), data = Microb_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8565 0.0515 74.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 8.029905)
##
## Null deviance: 528.4 on 63 degrees of freedom
## Residual deviance: 528.4 on 63 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Plot model assumption
plot_model(model3, type = "eff", terms = "SPECIES",
colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
theme(legend.text = element_text(face = "italic"))
Linear models with Poisson distribution (Phylogenetic q = 0)
With Microbial Diversity (q1)
divs$phyDiet<- divs$pq0
#Declare as factors
SPECIES <- as.factor(divs$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(divs$Season) # Dry and Rainy
SEX <- as.factor(divs$Sex) # Male and Female
#Standardize continous independent variables
SVL <- rescale(divs$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(divs$Elevation, binary.inputs = "center") # Taken as m a.s.l.
SEQDEPTH <- rescale(divs$SeqDepth, binary.inputs = "center")
## UPDATE FUNCTION (RUN ANOVA)
# 68 samples (diet / microbiota)
m1 <- glm(q1 ~ phyDiet*SPECIES+SEASON,
family = quasipoisson(link = "log"),
data = divs)
summary(m1)##
## Call:
## glm(formula = q1 ~ phyDiet * SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.661408 0.206567 22.566 < 2e-16 ***
## phyDiet -0.005446 0.069609 -0.078 0.93791
## SPECIESSceloporus bicanthalis -0.042498 0.236501 -0.180 0.85803
## SPECIESSceloporus grammicus -0.262152 0.243544 -1.076 0.28628
## SPECIESSceloporus spinosus -0.013347 0.247408 -0.054 0.95717
## SEASONRainy 0.276381 0.088922 3.108 0.00294 **
## phyDiet:SPECIESSceloporus bicanthalis 0.050503 0.084691 0.596 0.55332
## phyDiet:SPECIESSceloporus grammicus -0.004840 0.086901 -0.056 0.95577
## phyDiet:SPECIESSceloporus spinosus -0.014879 0.082267 -0.181 0.85711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.85698)
##
## Null deviance: 944.88 on 65 degrees of freedom
## Residual deviance: 619.92 on 57 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = q1 ~ phyDiet + SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.651537 0.133638 34.807 < 2e-16 ***
## phyDiet 0.001562 0.026733 0.058 0.95360
## SPECIESSceloporus bicanthalis 0.073228 0.106615 0.687 0.49483
## SPECIESSceloporus grammicus -0.274588 0.113016 -2.430 0.01812 *
## SPECIESSceloporus spinosus -0.062852 0.109585 -0.574 0.56842
## SEASONRainy 0.264348 0.086185 3.067 0.00324 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.49499)
##
## Null deviance: 944.88 on 65 degrees of freedom
## Residual deviance: 633.01 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.65694 0.09566 48.684 < 2e-16 ***
## SPECIESSceloporus bicanthalis 0.07275 0.10543 0.690 0.49280
## SPECIESSceloporus grammicus -0.27534 0.11136 -2.472 0.01622 *
## SPECIESSceloporus spinosus -0.06223 0.10818 -0.575 0.56720
## SEASONRainy 0.26231 0.07814 3.357 0.00136 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.32373)
##
## Null deviance: 944.88 on 65 degrees of freedom
## Residual deviance: 633.05 on 61 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.55181 0.06312 72.115 < 2e-16 ***
## SEASONRainy 0.32085 0.07995 4.013 0.00016 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 11.33081)
##
## Null deviance: 944.88 on 65 degrees of freedom
## Residual deviance: 758.28 on 64 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q1 ~ SPECIES, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.82632 0.08573 56.296 < 2e-16 ***
## SPECIESSceloporus bicanthalis 0.09084 0.11315 0.803 0.42514
## SPECIESSceloporus grammicus -0.33142 0.11834 -2.801 0.00679 **
## SPECIESSceloporus spinosus -0.08424 0.11604 -0.726 0.47057
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 11.91967)
##
## Null deviance: 944.88 on 65 degrees of freedom
## Residual deviance: 751.44 on 62 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Plot model assumption
plot_model(m3, type = "eff", terms = "SPECIES",
colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
theme(legend.text = element_text(face = "italic"))
plot_model(m3, type = "eff", terms = c("SPECIES", "SEASON"),
colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
theme_classic() +
theme(legend.text = element_text(face = "italic")) +
theme(axis.text.x = element_text(size = 11))
With Microbial Diversity (q2)
model1 <- glm(q2 ~ phyDiet*SPECIES+SEASON,
family = quasipoisson(link = "log"),
data = divs)
summary(model1)##
## Call:
## glm(formula = q2 ~ phyDiet * SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.878108 0.280118 13.845 <2e-16 ***
## phyDiet -0.042399 0.096811 -0.438 0.663
## SPECIESSceloporus bicanthalis 0.005418 0.319419 0.017 0.987
## SPECIESSceloporus grammicus -0.258721 0.327267 -0.791 0.432
## SPECIESSceloporus spinosus 0.045784 0.335106 0.137 0.892
## SEASONRainy 0.112807 0.115877 0.974 0.334
## phyDiet:SPECIESSceloporus bicanthalis 0.101063 0.114513 0.883 0.381
## phyDiet:SPECIESSceloporus grammicus 0.034825 0.117544 0.296 0.768
## phyDiet:SPECIESSceloporus spinosus -0.026253 0.113567 -0.231 0.818
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.702367)
##
## Null deviance: 568.65 on 65 degrees of freedom
## Residual deviance: 447.64 on 57 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = q2 ~ phyDiet + SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.82070 0.17848 21.407 <2e-16 ***
## phyDiet -0.01355 0.03571 -0.379 0.7058
## SPECIESSceloporus bicanthalis 0.24687 0.14407 1.714 0.0918 .
## SPECIESSceloporus grammicus -0.17655 0.15307 -1.153 0.2533
## SPECIESSceloporus spinosus -0.04872 0.15269 -0.319 0.7508
## SEASONRainy 0.08802 0.11442 0.769 0.4448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.670608)
##
## Null deviance: 568.65 on 65 degrees of freedom
## Residual deviance: 467.11 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SPECIES + SEASON, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.77434 0.12992 29.052 <2e-16 ***
## SPECIESSceloporus bicanthalis 0.24998 0.14286 1.750 0.0852 .
## SPECIESSceloporus grammicus -0.17020 0.15117 -1.126 0.2646
## SPECIESSceloporus spinosus -0.05479 0.15090 -0.363 0.7178
## SEASONRainy 0.10615 0.10332 1.027 0.3083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.57)
##
## Null deviance: 568.65 on 65 degrees of freedom
## Residual deviance: 468.22 on 61 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ SPECIES, family = quasipoisson(link = "log"),
## data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.84098 0.11157 34.427 <2e-16 ***
## SPECIESSceloporus bicanthalis 0.25751 0.14236 1.809 0.0753 .
## SPECIESSceloporus grammicus -0.19303 0.14920 -1.294 0.2006
## SPECIESSceloporus spinosus -0.06383 0.15031 -0.425 0.6725
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.536168)
##
## Null deviance: 568.65 on 65 degrees of freedom
## Residual deviance: 476.25 on 62 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = q2 ~ 1, family = quasipoisson(link = "log"), data = divs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.84306 0.05214 73.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 8.3741)
##
## Null deviance: 568.65 on 65 degrees of freedom
## Residual deviance: 568.65 on 65 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4# Plot model assumption
plot_model(model3, type = "eff", terms = "SPECIES",
colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
theme(legend.text = element_text(face = "italic"))
III. BETA DIVERSITY ANALYSES
Install hilldiv2
Loading libraries
Calculate beta diversity
## The estimated time for calculating the 4465 pairwise combinations is 5 seconds.## [1] "q1S" "q1C" "q1U" "q1V" "q2S" "q2C" "q2U" "q2V"## [1] 0.6500844 0.6843612 0.7557878 0.1858421 0.5226403 0.6284240- S (Jaccard-type turnover): it quantifies the normalised species turnover rate in a sample relative to the total pool of samples (i.e., gamma diversity).
- V (Sørensen-type turnover): it quantifies the normalised species turnover rate in a sample relative to one sample (i.e., alpha diversity).
- U (Jaccard-type overlap-complement): it quantifies the proportion of non-shared species in the total pool of samples. Therefore, this metric quantifies dissimilarity from the perspective of the pool of samples.
- C (Sørensen-type overlap-complement): it quantifies the effective average proportion of non-shared OTUs/ASVs/MAGs in samples. Therefore, this metric quantifies dissimilarity from the perspective of a single sample.
Taxonomic beta diversity at order q1
hill_pair_dis_nmds <- hill_pair_dis$q1S %>% #select the distance object based on dissimilarity metric S
metaMDS(.,trymax = 500, k = 2, verbosity = FALSE) %>%
vegan:::scores.metaMDS() %>%
as_tibble(., rownames = "sample")
hill_pair_dis_nmds <- hill_pair_dis_nmds %>%
left_join(Metadata, by = join_by(sample == SampleID)) %>%
group_by(Species, Season) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup()Plot ordination q1
NMDS_Plot_q1 <- ggplot(hill_pair_dis_nmds, aes(x = NMDS1, y = NMDS2,
color = Species, shape = Season)) +
geom_point(size = 2) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.2) +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
theme_classic() +
theme(legend.text = element_text(face = "italic")) +
scale_colour_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))
print(NMDS_Plot_q1)
Run PERMANONA
set.seed(123)
dist <- hill_pair_dis$q1S
perm_q1 <- adonis2(dist ~ Species*Season+Elevation+SVL+SeqDepth,
data = Metadata, permutations = 999)
print(perm_q1)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist ~ Species * Season + Elevation + SVL + SeqDepth, data = Metadata, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## Species 3 6.6636 0.22329 10.2239 0.001 ***
## Season 1 1.7738 0.05944 8.1645 0.001 ***
## Elevation 1 0.2071 0.00694 0.9531 0.566
## SVL 1 0.3386 0.01134 1.5583 0.053 .
## SeqDepth 1 0.2866 0.00960 1.3191 0.118
## Species:Season 3 2.3235 0.07786 3.5648 0.001 ***
## Residual 84 18.2495 0.61153
## Total 94 29.8426 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Run pairwise perMANOVA
Metadata$species_season <- paste(Metadata$Species, Metadata$Season, sep = "_")
pw_adonis_q1<-pairwise.adonis(dist,factors=Metadata$species_season)
intrasp_compar <- pw_adonis_q1[grep("^([^_]+)_Dry vs \\1_Rainy", pw_adonis_q1$pairs), ]
print(intrasp_compar)## pairs Df SumsOfSqs
## 5 Sceloporus grammicus_Dry vs Sceloporus grammicus_Rainy 1 1.2919702
## 10 Sceloporus spinosus_Dry vs Sceloporus spinosus_Rainy 1 0.9619991
## 17 Sceloporus aeneus_Dry vs Sceloporus aeneus_Rainy 1 0.9720931
## 22 Sceloporus bicanthalis_Dry vs Sceloporus bicanthalis_Rainy 1 0.9324159
## F.Model R2 p.value p.adjusted sig
## 5 5.919983 0.1854632 0.001 0.028 .
## 10 3.704441 0.1631593 0.001 0.028 .
## 17 4.424882 0.1811629 0.001 0.028 .
## 22 5.104701 0.1883327 0.001 0.028 .Plot of intraspecific distances
compar_intrasp<- c("Sceloporus grammicus_Dry_vs_Sceloporus grammicus_Rainy", "Sceloporus aeneus_Dry_vs_Sceloporus aeneus_Rainy", "Sceloporus bicanthalis_Dry_vs_Sceloporus bicanthalis_Rainy", "Sceloporus spinosus_Dry_vs_Sceloporus spinosus_Rainy")
dist.q1 <- 1-(hill_pair_dis$q1S) %>% as.matrix()
#convert matrix into df
dist.q1[lower.tri(dist.q1)] <- NA
dist.q1.df <- as.matrix(dist.q1) %>% melt(
varnames = c("site1", "site2"))%>% drop_na() %>% filter(!value==0)
#Join beta diversity vakues with metadata
dist.q1.meta<- dist.q1.df %>%
inner_join(Metadata, by = c("site1"="SampleID")) %>%
inner_join(Metadata, by = c("site2"="SampleID"))
#Subset to obtain specific comparisons
dist.q1.meta.filt <- dist.q1.meta %>% unite(
"compar_intraspp", c("species_season.x", "species_season.y"), sep="_vs_", remove=F) %>% filter(compar_intraspp %in% compar_intrasp)
# Plotting
#Setting the order of the species
dist.q1.meta.filt$Species.order<- factor(dist.q1.meta.filt$Species.x,
levels = c("Sceloporus aeneus", "Sceloporus bicanthalis", "Sceloporus grammicus", "Sceloporus spinosus"))
# Kruskal wallis test
set.seed(1234)
kruskal.t.q1 <- dist.q1.meta.filt %>% kruskal_test(value ~ Species.order)
Post_hoc.q1 <- dist.q1.meta.filt %>% dunn_test(value ~ Species.order, p.adjust.method = "bonferroni")
Post_hoc.q1 <- Post_hoc.q1 %>%
add_xy_position(x = "Species.order")
# Boxplot
q1_intrasp<-dist.q1.meta.filt %>%
ggboxplot(x = "Species.order", y="value", fill = "Species.order")+
ylab("Seasonal turnover of common ASVs (q=1)")+
xlab(element_blank())+
scale_fill_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))+
theme_classic() +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank()) +
theme(legend.text = element_text(face = "italic"))+
stat_pvalue_manual(Post_hoc.q1,
hide.ns = TRUE,
step.increase = 0.2) +
labs(
subtitle = get_test_label(kruskal.t.q1,
detailed = TRUE),
caption = get_pwc_label(Post_hoc.q1))+
scale_y_continuous(limits = c(0,1))
q1_intrasp
Run BETADISPER analysis
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.06544 0.0218120 4.1684 999 0.008 **
## Residuals 91 0.47618 0.0052328
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.06544 0.0218120 4.1684 999 0.013 *
## Residuals 91 0.47618 0.0052328
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Sceloporus aeneus Sceloporus bicanthalis
## Sceloporus aeneus 0.17700000
## Sceloporus bicanthalis 0.16097306
## Sceloporus grammicus 0.66564466 0.01087924
## Sceloporus spinosus 0.11780891 0.00011976
## Sceloporus grammicus Sceloporus spinosus
## Sceloporus aeneus 0.66800000 0.122
## Sceloporus bicanthalis 0.01700000 0.001
## Sceloporus grammicus 0.091
## Sceloporus spinosus 0.09456538
Taxonomic beta diversity at order q2
hill_pair_dis_nmds_q2 <- hill_pair_dis$q2S %>% #select the distance object based on dissimilarity metric S
metaMDS(.,trymax = 500, k = 2, verbosity = FALSE) %>%
vegan:::scores.metaMDS() %>%
as_tibble(., rownames = "sample")
hill_pair_dis_nmds_2 <- hill_pair_dis_nmds_q2 %>%
left_join(Metadata, by = join_by(sample == SampleID)) %>%
group_by(Species, Season) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup()Plot ordination q2
NMDS_Plot_q2 <- ggplot(hill_pair_dis_nmds_2, aes(x = NMDS1, y = NMDS2,
color = Species, shape = Season)) +
geom_point(size = 2) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.2) +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
theme_classic() +
theme(legend.text = element_text(face = "italic")) +
scale_colour_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))
print(NMDS_Plot_q2)
Run PERMANOVA
set.seed(123)
dist2 <- hill_pair_dis$q2S
perm_q2 <- adonis2(dist2 ~ Species*Season+Elevation+SVL+SeqDepth,
data = Metadata, permutations = 999)
print(perm_q2)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist2 ~ Species * Season + Elevation + SVL + SeqDepth, data = Metadata, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## Species 3 6.665 0.19570 8.7137 0.001 ***
## Season 1 2.129 0.06251 8.3503 0.001 ***
## Elevation 1 0.226 0.00664 0.8867 0.659
## SVL 1 0.393 0.01155 1.5431 0.067 .
## SeqDepth 1 0.294 0.00865 1.1548 0.270
## Species:Season 3 2.933 0.08612 3.8346 0.001 ***
## Residual 84 21.418 0.62884
## Total 94 34.060 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Run pairwise perMANOVA
pw_adonis_q2<-pairwise.adonis(dist2,factors=Metadata$species_season)
intrasp_compar_q2 <- pw_adonis_q2[grep("^([^_]+)_Dry vs \\1_Rainy", pw_adonis_q2$pairs), ]
print(intrasp_compar_q2)## pairs Df SumsOfSqs
## 5 Sceloporus grammicus_Dry vs Sceloporus grammicus_Rainy 1 1.3213474
## 10 Sceloporus spinosus_Dry vs Sceloporus spinosus_Rainy 1 0.9123456
## 17 Sceloporus aeneus_Dry vs Sceloporus aeneus_Rainy 1 1.2367957
## 22 Sceloporus bicanthalis_Dry vs Sceloporus bicanthalis_Rainy 1 1.6336642
## F.Model R2 p.value p.adjusted sig
## 5 5.294296 0.1691777 0.001 0.028 .
## 10 2.860417 0.1308492 0.001 0.028 .
## 17 4.641605 0.1883646 0.001 0.028 .
## 22 8.146037 0.2702192 0.001 0.028 .Plot of intraspecific distances
compar_intrasp<- c("Sceloporus grammicus_Dry_vs_Sceloporus grammicus_Rainy", "Sceloporus aeneus_Dry_vs_Sceloporus aeneus_Rainy", "Sceloporus bicanthalis_Dry_vs_Sceloporus bicanthalis_Rainy", "Sceloporus spinosus_Dry_vs_Sceloporus spinosus_Rainy")
dist.q2 <- 1-(hill_pair_dis$q2S) %>% as.matrix()
#convert matrix into df
dist.q2[lower.tri(dist.q2)] <- NA
dist.q2.df <- as.matrix(dist.q2) %>% melt(
varnames = c("site1", "site2"))%>% drop_na() %>% filter(!value==0)
#Join beta diversity vakues with metadata
dist.q2.meta<- dist.q2.df %>%
inner_join(Metadata, by = c("site1"="SampleID")) %>%
inner_join(Metadata, by = c("site2"="SampleID"))
#Subset to obtain specific comparisons
dist.q2.meta.filt <- dist.q2.meta %>% unite(
"compar_intraspp", c("species_season.x", "species_season.y"), sep="_vs_", remove=F) %>% filter(compar_intraspp %in% compar_intrasp)
# Plotting
#Setting the order of the species
dist.q2.meta.filt$Species.order<- factor(dist.q2.meta.filt$Species.x,
levels = c("Sceloporus aeneus", "Sceloporus bicanthalis", "Sceloporus grammicus", "Sceloporus spinosus"))
# Kruskal wallis test
set.seed(1234)
kruskal.t.q2 <- dist.q2.meta.filt %>% kruskal_test(value ~ Species.order)
Post_hoc.q2 <- dist.q2.meta.filt %>% dunn_test(value ~ Species.order, p.adjust.method = "bonferroni")
Post_hoc.q2 <- Post_hoc.q2 %>%
add_xy_position(x = "Species.order")
q2_intrasp<-dist.q2.meta.filt %>%
ggboxplot(x = "Species.order", y="value", fill = "Species.order")+
ylab("Seasonal turnover of dominant ASVs (q=2)")+
xlab(element_blank())+
scale_fill_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))+
theme_classic() +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
legend.title = element_blank()) +
theme(legend.text = element_text(face = "italic"))+
stat_pvalue_manual(Post_hoc.q2,
hide.ns = TRUE,
step.increase = 0.1) +
labs(
subtitle = get_test_label(kruskal.t.q2,
detailed = TRUE),
caption = get_pwc_label(Post_hoc.q2))+
scale_y_continuous(limits = c(0,1))
q2_intrasp
Run BETADISPER analysis
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.07338 0.0244593 2.5963 999 0.043 *
## Residuals 91 0.85729 0.0094208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.07338 0.0244593 2.5963 999 0.067 .
## Residuals 91 0.85729 0.0094208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Sceloporus aeneus Sceloporus bicanthalis
## Sceloporus aeneus 0.2130000
## Sceloporus bicanthalis 0.2201605
## Sceloporus grammicus 0.6839133 0.2468309
## Sceloporus spinosus 0.2294011 0.0059616
## Sceloporus grammicus Sceloporus spinosus
## Sceloporus aeneus 0.6960000 0.232
## Sceloporus bicanthalis 0.2560000 0.007
## Sceloporus grammicus 0.030
## Sceloporus spinosus 0.0297543
IV. MANTEL Y PROCRUSTES ANALYSIS
Load libraries
Loading files
This analysis is based on 68 samples obtained from both datasets (diet-microbiota)
metadata <- read.csv("../data/Procrustes_metadata.csv", header = TRUE, check.names = FALSE,
row.names = 1)
diet <- read.csv("../data/Procrustes_Diet.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota <- read.csv("../data/Procrustes_Microbiota.csv", header = TRUE,
check.names = FALSE, row.names = 1)
#microbiota <- microbiota %>% filter(rowSums(across(where(is.numeric)))!=0) %>% t()Calculate distances matrices
# Distance matrix of bacterial microbiota
jaccard_mic_1 <- hillpair(data = microbiota, q = 1) # distance matrix at order q1
microbiota_jaccard_1 <- as.dist(jaccard_mic_1$S)
jaccard_mic_2 <- hillpair(data = microbiota, q = 2) # distance matrix at order q1
microbiota_jaccard_2 <- as.dist(jaccard_mic_2$S)
# Distance matrix of diet
diet_jaccard <- vegdist(diet, method = "jaccard")Mantel test
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard, ydis = microbiota_jaccard_1, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.206
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0349 0.0473 0.0560 0.0687
## Permutation: free
## Number of permutations: 999##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard, ydis = microbiota_jaccard_2, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1948
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0387 0.0511 0.0610 0.0764
## Permutation: free
## Number of permutations: 999Make pcoas
# PCoA of microbiota q=1
microbiotaPCoA_1 <- as.data.frame(cmdscale(microbiota_jaccard_1))
plot(microbiotaPCoA_1)
# PCoA of microbiota q=2
microbiotaPCoA_2 <- as.data.frame(cmdscale(microbiota_jaccard_2))
plot(microbiotaPCoA_2)
Procrustes analysis of all species
Procrustes analysis of all species at order q=1
procrust <- procrustes(X = dietPCoA, Y = microbiotaPCoA_1, scale=TRUE,
symmetric = TRUE)
pro_test <- protest(dietPCoA, microbiotaPCoA_1, permutations = 9999)
pro_test##
## Call:
## protest(X = dietPCoA, Y = microbiotaPCoA_1, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.8268
## Correlation in a symmetric Procrustes rotation: 0.4161
## Significance: 1e-04
##
## Permutation: free
## Number of permutations: 9999
eigen <- sqrt(procrust$svd$d)
percent_var <- signif(eigen/sum(eigen), 4)*100
beta_pro <- data.frame(procrust$X)
head(beta_pro)beta_pro$UserName <- rownames(beta_pro)
beta_pro$type <- "Diet"
beta_pro$seasons<-metadata[,5]
beta_pro$species<-metadata[,1]
head(beta_pro)trans_pro$UsarName <- rownames(trans_pro)
trans_pro$type <- "Microbiota"
trans_pro$seasons=metadata[,5]
trans_pro$species=metadata[,1]
head(trans_pro)colnames(trans_pro) <- colnames(beta_pro)
toplot <- data.frame(rbind(beta_pro,trans_pro))
pval <- signif(pro_test$signif, 1)
#col1 = rgb(250/255,60/255,60/255)
#col4 = rgb(0/255,200/255,200/255)
#col8 = rgb(160/255,0/255,200/255)
#col10 = rgb(0/255,160/255,255/255)
# "#999999","#E69F00","#56B4E9", "#009E73","#F0E442Graph diet and microbiota
procrustes_new <- ggplot(toplot) +
geom_point(size = 1.8, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#999999","#E69F00","#56B4E9",
"#009E73")) +
theme_classic() +
scale_x_continuous(limits = c(-0.1, 0.3)) +
scale_y_continuous(limits = c(-0.2, 0.2)) +
geom_line(aes(x= V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
annotate("text", x = 0.2, y = -0.2,
label = paste0("p-value=", pval), size = 4) +
xlab(paste0("PCoA 1 (",percent_var[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var[2],"%)"))
print(procrustes_new)
Procrustes analysis of all species at order q=2
procrust <- procrustes(X = dietPCoA, Y = microbiotaPCoA_2, scale=TRUE,
symmetric = TRUE)
pro_test <- protest(dietPCoA, microbiotaPCoA_2, permutations = 9999)
pro_test##
## Call:
## protest(X = dietPCoA, Y = microbiotaPCoA_2, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.8177
## Correlation in a symmetric Procrustes rotation: 0.4269
## Significance: 1e-04
##
## Permutation: free
## Number of permutations: 9999
eigen <- sqrt(procrust$svd$d)
percent_var <- signif(eigen/sum(eigen), 4)*100
beta_pro <- data.frame(procrust$X)
head(beta_pro)beta_pro$UserName <- rownames(beta_pro)
beta_pro$type <- "Diet"
beta_pro$seasons<-metadata[,5]
beta_pro$species<-metadata[,1]
head(beta_pro)trans_pro$UsarName <- rownames(trans_pro)
trans_pro$type <- "Microbiota"
trans_pro$seasons=metadata[,5]
trans_pro$species=metadata[,1]
head(trans_pro)Graph diet and microbiota
procrustes_new <- ggplot(toplot) +
geom_point(size = 1.8, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#999999","#E69F00","#56B4E9",
"#009E73")) +
theme_classic() +
scale_x_continuous(limits = c(-0.1, 0.3)) +
scale_y_continuous(limits = c(-0.2, 0.2)) +
geom_line(aes(x= V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
annotate("text", x = 0.2, y = -0.2,
label = paste0("p-value=", pval), size = 4) +
xlab(paste0("PCoA 1 (",percent_var[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var[2],"%)"))
print(procrustes_new)
Procrustes analysis by lizard species
Procrsutes by lizard species were done at q=1 as similar patterns were observed in the general procrustes at q=2.
Procrustes of Sceloporus aeneus
#loading data
aeneus_metadata <- subset(metadata, Species_Scel == "Sceloporus aeneus")
diet_aeneus <- read.csv("../data/Proc_Diet_aeneus.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_aeneus <- read.csv("../data/Proc_Mic_aeneus.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_a <- hillpair(data = microbiota_aeneus, q = 1)
microbiota_jaccard_a <- as.dist(jaccard_mic_a$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_a <- vegdist(diet_aeneus, method = "jaccard")
diet_jaccard_a <- as.dist(diet_jaccard_a)
# Mantel test
set.seed(123)
mantel(diet_jaccard_a, microbiota_jaccard_a, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_a, ydis = microbiota_jaccard_a, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.3412
## Significance: 0.031
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.196 0.283 0.348 0.399
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_a <- procrustes(X = dietPCoA_a, Y = microbiotaPCoA_a, scale=TRUE,
symmetric = TRUE)
pro_test_a <- protest(dietPCoA_a, microbiotaPCoA_a, permutations = 9999)
pro_test_a##
## Call:
## protest(X = dietPCoA_a, Y = microbiotaPCoA_a, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.7374
## Correlation in a symmetric Procrustes rotation: 0.5125
## Significance: 0.0734
##
## Permutation: free
## Number of permutations: 9999
eigen_a <- sqrt(procrust_a$svd$d)
percent_var_a <- signif(eigen_a/sum(eigen_a), 4)*100
beta_pro_a <- data.frame(procrust_a$X)
head(beta_pro_a)beta_pro_a$UserName <- rownames(beta_pro_a)
beta_pro_a$type <- "Diet"
beta_pro_a$seasons <- aeneus_metadata[,5]
beta_pro_a$species <- aeneus_metadata[,1]
head(beta_pro_a)trans_pro_a$UsarName <- rownames(trans_pro_a)
trans_pro_a$type <- "Microbiota"
trans_pro_a$seasons <- aeneus_metadata[,5]
trans_pro_a$species <- aeneus_metadata[,1]
head(trans_pro_a)colnames(trans_pro_a) <- colnames(beta_pro_a)
toplot_a <- data.frame(rbind(beta_pro_a, trans_pro_a))
pval_a <- signif(pro_test_a$signif, 1)
# Graph diet and microbiota
procreustes_aeneus_new <- ggplot(toplot_a) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#999999")) +
theme_classic() +
scale_x_continuous(limits = c(-0.2, 0.4)) +
scale_y_continuous(limits = c(-0.2, 0.4)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.13, label = paste0("p-value=", pval_a),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_a[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_a[2],"%)"))
print(procreustes_aeneus_new)
Procrustes of Sceloporus bicanthalis
#Loading data
bica_metadata <- subset(metadata, Species_Scel == "Sceloporus bicanthalis")
diet_bica <- read.csv("../data/Proc_Diet_bica.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_bica <- read.csv("../data/Proc_Mic_bica.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_b <- hillpair(data = microbiota_bica, q = 1)
microbiota_jaccard_b <- as.dist(jaccard_mic_b$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_b <- vegdist(diet_bica, method = "jaccard")
diet_jaccard_b <- as.dist(diet_jaccard_b)
# Mantel test
set.seed(123)
mantel(diet_jaccard_b, microbiota_jaccard_b, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_b, ydis = microbiota_jaccard_b, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.3353
## Significance: 0.006
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.135 0.182 0.219 0.284
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_b <- procrustes(X = dietPCoA_b, Y = microbiotaPCoA_b, scale=TRUE,
symmetric = TRUE)
pro_test_b <- protest(dietPCoA_b, microbiotaPCoA_b, permutations = 9999)
pro_test_b##
## Call:
## protest(X = dietPCoA_b, Y = microbiotaPCoA_b, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.7286
## Correlation in a symmetric Procrustes rotation: 0.521
## Significance: 0.0226
##
## Permutation: free
## Number of permutations: 9999
eigen_b <- sqrt(procrust_b$svd$d)
percent_var_b <- signif(eigen_b/sum(eigen_b), 4)*100
beta_pro_b <- data.frame(procrust_b$X)
head(beta_pro_b)beta_pro_b$UserName <- rownames(beta_pro_b)
beta_pro_b$type <- "Diet"
beta_pro_b$seasons <- bica_metadata[,5]
beta_pro_b$species <- bica_metadata[,1]
head(beta_pro_b)trans_pro_b$UsarName <- rownames(trans_pro_b)
trans_pro_b$type <- "Microbiota"
trans_pro_b$seasons <- bica_metadata[,5]
trans_pro_b$species <- bica_metadata[,1]
head(trans_pro_b)colnames(trans_pro_b) <- colnames(beta_pro_b)
toplot_b <- data.frame(rbind(beta_pro_b, trans_pro_b))
pval_b <- signif(pro_test_b$signif, 1)
# graph diet and microbiota
procreustes_bica_new <- ggplot(toplot_b) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#E69F00")) +
theme_classic() +
#scale_x_continuous(limits = c(-0.25, 0.4)) +
#scale_y_continuous(limits = c(-0.2, 0.25)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_b),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_b[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_b[2],"%)"))
print(procreustes_bica_new)
Procrustes: Sceloporus grammicus
# Loading data
gram_metadata <- subset(metadata, Species_Scel == "Sceloporus grammicus")
diet_gram <- read.csv("../data/Proc_Diet_gram.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_gram <- read.csv("../data/Proc_Mic_gram.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_g <- hillpair(data = microbiota_gram, q = 1)
microbiota_jaccard_g <- as.dist(jaccard_mic_g$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_g <- vegdist(diet_gram, method = "jaccard")
diet_jaccard_g <- as.dist(diet_jaccard_g)
# Mantel test
set.seed(123)
mantel(diet_jaccard_g, microbiota_jaccard_g, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_g, ydis = microbiota_jaccard_g, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.4457
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0905 0.1235 0.1568 0.1910
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_g <- procrustes(X = dietPCoA_g, Y = microbiotaPCoA_g, scale=TRUE,
symmetric = TRUE)
pro_test_g <- protest(dietPCoA_g, microbiotaPCoA_g, permutations = 9999)
pro_test_g##
## Call:
## protest(X = dietPCoA_g, Y = microbiotaPCoA_g, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.5908
## Correlation in a symmetric Procrustes rotation: 0.6397
## Significance: 1e-04
##
## Permutation: free
## Number of permutations: 9999
eigen_g <- sqrt(procrust_g$svd$d)
percent_var_g <- signif(eigen_g/sum(eigen_g), 4)*100
beta_pro_g <- data.frame(procrust_g$X)
head(beta_pro_g)beta_pro_g$UserName <- rownames(beta_pro_g)
beta_pro_g$type <- "Diet"
beta_pro_g$seasons <- gram_metadata[,5]
beta_pro_g$species <- gram_metadata[,1]
head(beta_pro_g)trans_pro_g$UsarName <- rownames(trans_pro_g)
trans_pro_g$type <- "Microbiota"
trans_pro_g$seasons <- gram_metadata[,5]
trans_pro_g$species <- gram_metadata[,1]
head(trans_pro_g)colnames(trans_pro_g) <- colnames(beta_pro_g)
toplot_g <- data.frame(rbind(beta_pro_g, trans_pro_g))
pval_g <- signif(pro_test_g$signif, 1)
# Graph diet and microbiota
procreustes_gram_new <- ggplot(toplot_g) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#56B4E9")) +
theme_classic() +
#scale_x_continuous(limits = c(-0.25, 0.4)) +
#scale_y_continuous(limits = c(-0.2, 0.25)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_b),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_g[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_g[2],"%)"))
print(procreustes_gram_new)
Procrustes: Sceloporus spinosus
spi_metadata <- subset(metadata, Species_Scel == "Sceloporus spinosus")
diet_spi <- read.csv("../data/Proc_Diet_spi.csv", header = TRUE, check.names = FALSE,
row.names = 1) %>% t()
microbiota_spi <- read.csv("../data/Proc_Mic_spi.csv", header = TRUE,
check.names = FALSE, row.names = 1)
# Distance matrix
jaccard_mic_s <- hillpair(data = microbiota_spi, q = 1)
microbiota_jaccard_s <- as.dist(jaccard_mic_s$S)
# Absence/presence data frame - Jaccard distance
diet_jaccard_s <- vegdist(diet_spi, method = "jaccard")
diet_jaccard_s <- as.dist(diet_jaccard_s)
# Mantel test
set.seed(123)
mantel(diet_jaccard_s, microbiota_jaccard_s, method = "pearson",
permutations = 999)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = diet_jaccard_s, ydis = microbiota_jaccard_s, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1073
## Significance: 0.151
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.141 0.190 0.224 0.267
## Permutation: free
## Number of permutations: 999
# Procrustes analysis
procrust_s <- procrustes(X = dietPCoA_s, Y = microbiotaPCoA_s, scale=TRUE,
symmetric = TRUE)
pro_test_s <- protest(dietPCoA_s, microbiotaPCoA_s, permutations = 9999)
pro_test_s##
## Call:
## protest(X = dietPCoA_s, Y = microbiotaPCoA_s, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.797
## Correlation in a symmetric Procrustes rotation: 0.4505
## Significance: 0.0789
##
## Permutation: free
## Number of permutations: 9999
eigen_s <- sqrt(procrust_s$svd$d)
percent_var_s <- signif(eigen_s/sum(eigen_s), 4)*100
beta_pro_s <- data.frame(procrust_s$X)
head(beta_pro_s)beta_pro_s$UserName <- rownames(beta_pro_s)
beta_pro_s$type <- "Diet"
beta_pro_s$seasons <- spi_metadata[,5]
beta_pro_s$species <- spi_metadata[,1]
head(beta_pro_s)trans_pro_s$UsarName <- rownames(trans_pro_s)
trans_pro_s$type <- "Microbiota"
trans_pro_s$seasons <- spi_metadata[,5]
trans_pro_s$species <- spi_metadata[,1]
head(trans_pro_s)colnames(trans_pro_s) <- colnames(beta_pro_s)
toplot_s <- data.frame(rbind(beta_pro_s, trans_pro_s))
pval_s <- signif(pro_test_s$signif, 1)
# Graph diet and microbiota
procreustes_spi_new <- ggplot(toplot_s) +
geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species,
shape = type)) +
scale_color_manual(values = c("#40A175")) +
theme_classic() +
#scale_x_continuous(limits = c(-0.25, 0.4)) +
#scale_y_continuous(limits = c(-0.2, 0.25)) +
geom_line(aes(x = V1, y = V2, group = UserName),
col = "darkgrey", alpha = 0.6,
size = 0.2) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10, colour = "black"),
legend.position = "right",
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
aspect.ratio = 1) +
theme(legend.text = element_text(face = "italic")) +
guides(color = guide_legend(ncol = 1)) +
#annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_s),
# size = 4) +
xlab(paste0("PCoA 1 (",percent_var_s[1],"%)")) +
ylab(paste0("PCoA 2 (",percent_var_s[2],"%)"))
print(procreustes_spi_new)
Procustes and mantel analysis considering phylogenetic comosition of diet
tree.diet <- read.tree("../data/tree.nwk")
comm<-(diet)
IDs= colnames(comm) # para obtener los nombres de las sp
tree.diet= ape::rtree(n=ncol(diet), tip.label = paste0(IDs)) # tener el árbol enraizado, con los nombres de los sitios y los nombres de las sp.
hill_phylo_parti_pairwise(comm,tree.diet, q = 0, show_warning = F) # Loading phylogenetic tree
phy_tree <- read_tree("../data/tree.nwk")
# Creating phyloseq object
SAM <- sample_data(metadata)
OTU <- otu_table(diet, taxa_are_rows=F)
phylo_physeq = phyloseq(OTU, SAM, phy_tree)
#Calculating phylogenetic distance
diet_dist_phy = distance(phylo_physeq, method = "uunifrac") %>% as.dist